# Python simulation to illustrate expectation, p-th moments, and sample vs theoretical behavior.
# This code will:
# 1. Simulate samples from Normal(0,1), Exponential(1), and Cauchy(0,1).
# 2. Compute cumulative (running) sample means and running p-th absolute moments (p=1,2,3).
# 3. Show how sample quantities behave as sample size grows (Law of Large Numbers / failure cases).
# 4. Present a small table of values at several sample sizes and final statistics (n = 100000).
#
# NOTE: plots use matplotlib (no seaborn) and each chart is a separate figure.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib
matplotlib.use(backend="TkAgg")

np.random.seed(2025)

max_n = 100_000
sizes_to_show = np.array([10, 100, 1_000, 5_000, 10_000, 20_000, 50_000, 100_000])

# generate samples
samples = {
    "Normal(0,1)": np.random.normal(loc=0.0, scale=1.0, size=max_n),
    "Exponential(1)": np.random.exponential(scale=1.0, size=max_n),  # mean = 1
    "Cauchy(0,1)": np.random.standard_cauchy(size=max_n)              # heavy-tailed, mean undefined
}

# We'll store running statistics and summary tables
running_means = {}
running_p_roots = {}  # keys: dist -> {p: running p-root array}
summary_at_sizes = {}
final_summary = []

p_list = [1, 2, 3]

for name, data in samples.items():
    # cumulative mean
    cumsum = np.cumsum(data)
    cummean = cumsum / np.arange(1, max_n + 1)
    running_means[name] = cummean

    # running p-th absolute moments root: (mean(|X|^p))^(1/p)
    running_p_roots[name] = {}
    for p in p_list:
        cum_abs_p = np.cumsum(np.abs(data) ** p)
        running_abs_p_mean = cum_abs_p / np.arange(1, max_n + 1)
        running_p_roots[name][p] = running_abs_p_mean ** (1.0 / p)

    # summary values at selected sizes
    rows = []
    for s in sizes_to_show:
        row = {"n": int(s), "running_mean": float(cummean[int(s) - 1])}
        for p in p_list:
            row[f"root_p{p}"] = float(running_p_roots[name][p][int(s) - 1])
        rows.append(row)
    summary_at_sizes[name] = pd.DataFrame(rows).set_index("n")

    # final overall statistics (n = max_n)
    final_summary.append({
        "dist": name,
        "n": max_n,
        "sample_mean": float(cummean[-1]),
        "root_p1": float(running_p_roots[name][1][-1]),
        "root_p2": float(running_p_roots[name][2][-1]),
        "root_p3": float(running_p_roots[name][3][-1])
    })

# Display tables for the user
print("=== 样本在不同 n 下的运行统计（部分 n） ===\n")
for name, df in summary_at_sizes.items():
    print(f"--- {name} ---")
    display(df)  # pandas DataFrame will be displayed nicely

final_df = pd.DataFrame(final_summary).set_index("dist")
print("\n=== n = 100000 时的最终统计 ===")
display(final_df)

# Now plot running sample mean for each distribution (separate figures)
for name, cummean in running_means.items():
    plt.figure(figsize=(8, 3.5))
    plt.semilogx(np.arange(1, max_n + 1), cummean)
    plt.axhline(0, linestyle='--', linewidth=0.6) if name.startswith("Normal") else None
    plt.xlabel("样本量 n (对数刻度)")
    plt.ylabel("运行样本均值 (1/n Σ X_i)")
    plt.title(f"样本均值随 n 增大变化 — {name}")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Plot running p-root (只画 p=1 与 p=2) for each distribution
for name, p_dict in running_p_roots.items():
    plt.figure(figsize=(8, 3.5))
    plt.semilogx(np.arange(1, max_n + 1), p_dict[1], label="root p=1 (mean absolute)")
    plt.semilogx(np.arange(1, max_n + 1), p_dict[2], label="root p=2 (RMS)")
    plt.xlabel("样本量 n (对数刻度)")
    plt.ylabel("运行 p 阶均值 根值")
    plt.title(f"p 阶均值根随 n 增大变化 — {name}")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# Extra: quick histograms to visualize tail behavior (first 20000 samples)
for name, data in samples.items():
    plt.figure(figsize=(6,3.2))
    subset = data[:20000]
    # clip extreme values for visualization (but don't change underlying stats)
    clipped = np.clip(subset, -50, 50)
    plt.hist(clipped, bins=200, density=True)
    plt.xlabel("样本值 (极端值被裁剪到[-50,50]以便可视化)")
    plt.ylabel("相对频率")
    plt.title(f"样本直方图 (前20000, 裁剪) — {name}")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

print("（说明：Cauchy 的理论均值和二阶矩均未定义，因此对应运行统计会表现出非常大的波动，样本均值不会稳定收敛。）")
